% 光源产生需要考虑以下几个问题：
%     光源的光强分布：高斯分布（与位置相关）
%     光源的波长分布：高斯分布（与位置无关）
%     考虑光源是一个时稳光源，即光强分布随时间呈均匀变化，则光源只需要产生光子在空间上的密度分布和频率分布。

% 考虑到螺旋相位板对光源的转换率。
% 首先每个位置点的Direction只与散射角和光源宽度相关。
% 光强是统计信息，需要由概率分布和电场叠加求得。每个点的波长，高斯或OAM属性，相位则由都概率分布获得。
% 光强统计需要基于电场叠加，其余的信息经过光子的模拟。
% 每个点的电场信息是高斯电场和OAM电场的叠加，高斯电场就是OAM电场在l=0时的特殊情况。两个电场的权重是由粒子进行转换的平均比例决定的。
% 有了每个点的电场信息，每个点的光强和相位都会知道。
% 波长的概率密度分布仍然不变，

% api:
%     getRadius：得到光束半径
%     getLambda：得到光子波长
%     getElectricField：得到电场信息
%     getLightIntensity：得到光强信息
%     getgetDirection：得到方向余弦
%     getPhase：得到相位信息
%     getInfo：得到光子所有信息，包含光子波长，方向余弦，相位信息（使用电场表示） 

% 前提假设：
%     波长随空间分布均匀，转换率随空间分布均匀。（即两个分布都与空间无关），假设波长是正态分布，转换率是常量。

classdef lightSource
    properties
        % 与光子的相干高斯分布相关的参数
        miu;               % 波长分布的均值
        half_length;       % 波长分布半波带长度
        sigma;             % 波长分布的标准差
        % 高斯光到OAM转化率
        trans_rate1;       % 透镜1转换率
        trans_rate2;       % 透镜2转化率
        % 与衍射相关参数
        diffraction_angle; % 衍射角 / rad
        % OAM特征参数
        l;                 % OAM l
        p;                 % OAM p
        % 光源尺寸特征
        r;                 % 光束的束腰半径 
        % 衍生量
        k;                 % 平均波数
        z;                 % 传输距离
        z_R;               % 平均瑞利长度
        r_Z;               % 光束在z位置的半径
    end

    methods(Static)
        % 拉盖尔多项式
        function r = laguerre(p, l, x)
            if p == 0
                r = 1;
            elseif p == 1
                r = 1+abs(l)-x;
            else
                r = (1/p)*((2*p+l-1-x).*laguerre(p-1,abs(l),x)-(p+l-1)*laguerre(p-2,abs(l),x));
            end
        end
    end

    methods
        % 构造函数
        function obj = lightSource()
            obj.miu = 275e-9;                                   % 275nm
            obj.half_length = 6e-9;                             % 6nm
            obj.sigma = sqrt((obj.half_length/2)^2/(2*log(2))); % 标准差
            % 转化相关参数
            obj.trans_rate1 = 0.9;                              % 第一个透镜的转化率
            obj.trans_rate2 = 0.9;                              % 第二个透镜的转化率
            % 与衍射相关参数
            obj.diffraction_angle = 0;%0.001;%0.0262;                              % 26.2 mrad
            % OAM特征参数
            obj.l = 3;                                          % OAM l
            obj.p = 0;                                          % OAM p
            % 光源尺寸特征
            obj.r = 3e-3;                                       % 3mm
            % 衍生量
            obj.k = 2*pi/obj.miu;                               % 平均波数
            obj.z = 0;                                          % 传输距离，光源处
            obj.z_R = pi*(obj.r^2)/obj.miu;                     % 平均瑞利长度
            obj.r_Z = obj.r*sqrt(1+(obj.z/obj.z_R)^2);          % 光束在z位置的半径
        end

        % 获取光束半径
        function r = getRadius(obj)
            r = obj.r;
        end

        % 获取只经过透镜1半径光强比例
        function res = getTransOnceScale(obj, rdu, NNtr)
            R = obj.r*3;
            [x,y] = meshgrid(linspace(-R, R, NNtr));
            [theta, radius] = cart2pol(x,y);
            % scale1 OAM+ (1-scale1) 高斯
            E1 = sqrt(2*factorial(0)/pi/(0+factorial(abs(0))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(0)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(0,abs(0),2*radius.^2/obj.r_Z^2).*exp(-1i*0*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*0+abs(0)+1)*atan(obj.z/obj.z_R));
            E2 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E = (1-obj.trans_rate1)*E1 + obj.trans_rate1*E2;
            I = E.*conj(E);  
            I = I/max(max(I));
            I_target = 0;
            for i = 0 : NNtr-1
                for j = 0 : NNtr-1
                    xx = -R + 2*R*i/(NNtr-1);
                    yy = -R + 2*R*j/(NNtr-1);
                    if xx*xx + yy*yy <= rdu*rdu
                        % 在半径内
                        I_target = I_target + I(i+1, j+1);
                    end
                end
            end
            res = I_target/sum(sum(I));
        end

        % 获取经过两个透镜半径光强比例
        function res = getTransTwiceScale(obj, rdu, NNtr)
            R = obj.r*3;
            [x,y] = meshgrid(linspace(-R, R, NNtr));
            [theta, radius] = cart2pol(x,y);
            % scale1*scale2 + (1-scale1)(1-scale2) 高斯；scale1(1-scale2) OAM+ (1-scale1)scale2 OAM-
            E1 = sqrt(2*factorial(0)/pi/(0+factorial(abs(0))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(0)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(0,abs(0),2*radius.^2/obj.r_Z^2).*exp(-1i*0*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*0+abs(0)+1)*atan(obj.z/obj.z_R));
            E2 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E3 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*-obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E = (obj.trans_rate1*obj.trans_rate2 + (1-obj.trans_rate1)*(1-obj.trans_rate2))*E1 + obj.trans_rate1*(1-obj.trans_rate2)*E2 + (1-obj.trans_rate1)*obj.trans_rate2*E3;
            I = E.*conj(E);
            I = I/max(max(I));
            I_target = 0;
            for i = 0 : NNtr-1
                for j = 0 : NNtr-1
                    xx = -R + 2*R*i/(NNtr-1);
                    yy = -R + 2*R*j/(NNtr-1);
                    if xx*xx + yy*yy <= rdu*rdu
                        % 在半径内
                        I_target = I_target + I(i+1, j+1);
                    end
                end
            end
            res = I_target/sum(sum(I));
        end

        % 获取坐标点的波长
        function lambda = getLambda(obj)
            % 有限制条件，
            % Box–Muller算法
            a = rand;     
            b = rand;     
            X = sqrt((-2)*log(a))*cos(2*pi*b);
            lambda = X*obj.sigma + obj.miu;
        end

        % 获取某点的电场
        function E = getElectricField(obj, theta, radius)
            % E1为未转换的高斯电场，E2为转换后的OAM电场，两者的比例由转换的概率分布决定。
            E1 = sqrt(2*factorial(0)/pi/(0+factorial(abs(0))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(0)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(0,abs(0),2*radius.^2/obj.r_Z^2).*exp(-1i*0*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*0+abs(0)+1)*atan(obj.z/obj.z_R));
            E2 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E3 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*-obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E = ((obj.trans_rate1*obj.trans_rate2) + (1-obj.trans_rate1)*(1-obj.trans_rate2))*E1 + obj.trans_rate1*(1-obj.trans_rate2)*E2 + (1-obj.trans_rate1)*obj.trans_rate2*E3;
        end
        
        % 获取坐标点的光强
        function i = getLightIntensity(obj, theta, radius)
            e = obj.getElectricField(theta, radius);
            i = e.*conj(e);
        end

        % 获取坐标点的方向余弦
        function [u_x, u_y, u_z] = getDirection(obj, theta, radius)
            if obj.diffraction_angle == 0
                u_x = 0;
                u_y = 0;
                u_z = 1;
            else
                % 得到隐含的发射原点的坐标
                ori_x = 0;
                ori_y = 0;
                ori_z = -obj.r/tan(obj.diffraction_angle);
                % 极坐标转x,y,z坐标
                X = radius*cos(theta);
                Y = radius*sin(theta);
                Z = 0;
                % 得到发射向量
                vec_x = X - ori_x;
                vec_y = Y - ori_y;
                vec_z = Z - ori_z;
                length = sqrt(vec_x*vec_x + vec_y*vec_y + vec_z*vec_z);
                u_x = vec_x/length;
                u_y = vec_y/length;
                u_z = vec_z/length;
            end
        end

        % 获取坐标点的相位信息(使用电场表示)
        function [m, n, s, E1, E2] = getPhase(obj, theta, radius)
            % 以trans_rate决定此光子属于高斯光子还是OAM光子，如果是高斯光子，相位为0，如果是OAM光子，则由电场决定相位初始值。
            % 
            X = rand;
            if X < (obj.trans_rate1*obj.trans_rate2) + (1-obj.trans_rate1)*(1-obj.trans_rate2) % 高斯光，相位初始为0
                E1 = 1;
                E2 = 0;
                m = [1,0,0];
                n = [0,1,0];
                s = [0,0,1];
            elseif X < (obj.trans_rate1*obj.trans_rate2) + (1-obj.trans_rate1)*(1-obj.trans_rate2) + obj.trans_rate1*(1-obj.trans_rate2) % 被转化，OAM+
                e = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
                    .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
                    .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
                fi = angle(e);
                ele = cos(fi) + sin(fi)*1i;
                m = [1,0,0];
                n = [0,1,0];
                s = [0,0,1];
                E1 = ele;
                E2 = 0;
            else % OAM-
                e = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
                    .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*-obj.l*theta).*exp(-1i*obj.k*obj.z)...
                    .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
                fi = angle(e);
                ele = cos(fi) + sin(fi)*1i;
                m = [1,0,0];
                n = [0,1,0];
                s = [0,0,1];
                E1 = ele;
                E2 = 0;
            end
        end

        % 获取坐标点光子信息
        function [lambda, u_x, u_y, u_z, m, n, s, E1, E2] = getInfo(obj, theta, radius)
            lambda = obj.getLambda();
            [u_x, u_y, u_z] = obj.getDirection(theta, radius);
            [m, n, s, E1, E2] = obj.getPhase(theta, radius);
        end

        % 获取光强图
        function obj = showLightIntensityDiagram(obj, R, NNtr)
            [x,y] = meshgrid(linspace(-R, R, NNtr));
            [theta, radius] = cart2pol(x,y);
            E1 = sqrt(2*factorial(0)/pi/(0+factorial(abs(0))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(0)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(0,abs(0),2*radius.^2/obj.r_Z^2).*exp(-1i*0*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*0+abs(0)+1)*atan(obj.z/obj.z_R));
            E2 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E3 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*-obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E = ((obj.trans_rate1*obj.trans_rate2) + (1-obj.trans_rate1)*(1-obj.trans_rate2))*E1 + obj.trans_rate1*(1-obj.trans_rate2)*E2 + (1-obj.trans_rate1)*obj.trans_rate2*E3;
            I = E.*conj(E);  
            I = I/max(max(I));
            figure;
            mesh(x*1e3, y*1e3, I);
            colormap(jet);
            colorbar;
            view(2);
            set(gca,'fontname','times new roman','fontsize',14);
            title(['{\itl}=+', num2str(obj.l), ',{\itp}=', num2str(obj.p), ',{\itz}=', num2str(obj.z), 'm处的光强'], 'fontname', '华文中宋', 'fontsize', 14);
            xlabel('\itx/mm', 'fontname', 'times new roman', 'fontsize', 14);
            ylabel('\ity/mm', 'fontname', 'times new roman', 'fontsize', 14);
            zlabel('归一化强度', 'fontname', '华文中宋', 'fontsize', 14);
        end

        % 获取相位图
        function obj = showPhaseDiagram(obj, R, NNtr)
            [x,y] = meshgrid(linspace(-R, R, NNtr));
            [theta, radius] = cart2pol(x,y);
            E1 = sqrt(2*factorial(0)/pi/(0+factorial(abs(0))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(0)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(0,abs(0),2*radius.^2/obj.r_Z^2).*exp(-1i*0*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*0+abs(0)+1)*atan(obj.z/obj.z_R));
            E2 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E3 = sqrt(2*factorial(obj.p)/pi/(obj.p+factorial(abs(obj.l))))*(1/obj.r_Z)*(sqrt(2)*radius/obj.r_Z).^abs(obj.l)...
            .*exp(-radius.^2/obj.r_Z^2).*obj.laguerre(obj.p,abs(obj.l),2*radius.^2/obj.r_Z^2).*exp(-1i*-obj.l*theta).*exp(-1i*obj.k*obj.z)...
            .*exp(-1i*obj.k*radius.^2*obj.z/2/(obj.z^2+obj.z_R^2))*exp(-1i*(2*obj.p+abs(obj.l)+1)*atan(obj.z/obj.z_R));
            E = ((obj.trans_rate1*obj.trans_rate2) + (1-obj.trans_rate1)*(1-obj.trans_rate2))*E1 + obj.trans_rate1*(1-obj.trans_rate2)*E2 + (1-obj.trans_rate1)*obj.trans_rate2*E3;
            phase = angle(E);   
            phase = 2*pi*phase/max(max(phase));
            figure;
            mesh(x*1e3, y*1e3, phase);
            colormap(gray);
            colorbar;
            view(2);
            set(gca,'fontname','times new roman','fontsize',14);
            title(['{\itl}=+',num2str(obj.l), ',{\itp}=',num2str(obj.p), ',{\itz}=',num2str(obj.z), 'm处的相位'], 'fontname', '华文中宋', 'fontsize', 14);
            xlabel('\itx/mm', 'fontname', 'times new roman', 'fontsize', 14);
            ylabel('\ity/mm', 'fontname', 'times new roman', 'fontsize', 14);
            zlabel('归一化强度', 'fontname', '华文中宋', 'fontsize', 14);
        end
    end
end